Hierarchical Modelling - Practical 3

Hierarchical modelling using JAGS

We will continue with our hemlock light and growth example. Recall from the slides that we are modelling the relationship between plant growth rate and light as a non-linear curve defined by the Michaelis-Menten function - our process model. The equation here also indexes [s] for site which we come to later in the practical:

\[a[s] * (x[i,s]-c) / ((a[s]/b)+(x[i,s]-c))\]

# hemlock <- read.csv("Hemlock-light-data-simple.csv")
hemlock <- read.csv(url("https://raw.githubusercontent.com/NERC-CEH/beem_data/main/Hemlock-light-data-simple.csv"))
str(hemlock)

'data.frame':   77 obs. of  2 variables:
 $ Light               : num  98.5 98.5 98.5 23.5 11.8 11.8 19.4 19.4 98.7 98.7 ...
 $ Observed.growth.rate: num  18.2 26 25.7 6.2 10 11.2 7.3 15.7 44.2 39.7 ...

head(hemlock)

  Light Observed.growth.rate
1  98.5                 18.2
2  98.5                 26.0
3  98.5                 25.7
4  23.5                  6.2
5  11.8                 10.0
6  11.8                 11.2

Let’s read in the data, and take a quick look at how it’s structured and also plot it. We take a rudimentary approach to modelling the curve first, using non-linear least squares regression (a non-linear version of the familiar OLS).

par(mfrow=c(1,1))
plot(hemlock$Light, hemlock$Observed.growth.rate)

# Quick and dirty non-linear least squares fit to visualise
attach(hemlock)
start<-list(a=40, b=2, c=6)
m1<-nls(Observed.growth.rate~a*(Light-c)/((a/b)+(Light-c)), start=start)
summary(m1)

Formula: Observed.growth.rate ~ a * (Light - c)/((a/b) + (Light - c))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a   38.493      4.364   8.820 3.64e-13 ***
b    1.735      0.498   3.485 0.000831 ***
c    4.740      2.067   2.293 0.024664 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.409 on 74 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 5.505e-06

a<-coef(m1)["a"]
b<-coef(m1)["b"]
c<-coef(m1)["c"]
plot(Light,Observed.growth.rate)
x<-Light
curve(a*(x-c)/((a/b)+(x-c)), add=TRUE, col="blue")

detach(hemlock)

Now lets move ahead with the Bayesian analysis for a single non-linear fit using the software JAGS and the R package rjags, which ‘speaks to’ JAGS via our model file. First we need to set up the conditions for the model - how we want to run the MCMC chains, the data etc. Then we call to JAGS using our model file (open it and take a closer look).

Something worth noting is that we have tended to refer to our variance (sigma) in our model specifications, but JAGS uses precision, the inverse of variance (which typically takes the notation tau).

jags_model <- "
# JAGS model

# y is the observed growth rate
# x is the measurement of light, L

# a is alpha, the max. growth rate at infinite light
# c is light level at which growth is zero
# b is rate at which curve tails off (gamma)

########################################################################################
model{
  ### likelihood
  # Data model
  for (i in 1:n)
  {
    y[i] ~ dnorm(mu[i], tau.o) #Note JAGS uses precision, not variance
  }  
  
  # process model
  for (i in 1:n)
  {
    mu[i] ~ dnorm(mu2[i], tau.p) #Note JAGS uses precision, not variance
    
    mu2[i] <- a * (x[i]-c) / ((a/b)+(x[i]-c))
    
  }
  
  # priors
  a ~ dgamma(0.01,0.01)
  c ~ dunif(-10,10)
  b ~ dgamma(0.01,0.01)
  
  sigma.o ~ dnorm(5, 1/(0.5*0.5)) ## assume prior knowledge of observation error (5 with SD of 0.5)
  sigma.p ~ dunif(0, 50)
  
  # derived quantities: convert precisions to standard deviations
  tau.o <- 1/(sigma.o * sigma.o)
  tau.p <- 1/(sigma.p * sigma.p)
  
} # end of model
#########################################################################################
"
writeLines(jags_model, con="jags_model.txt")

library(rjags)

Warning: package 'rjags' was built under R version 4.4.3

Loading required package: coda

Warning: package 'coda' was built under R version 4.4.3

Linked to JAGS 4.3.1

Loaded modules: basemod,bugs

# Initialize the MCMC chains. This needs to be a list, so we define a function that calls the list. Anything that is to be ESTIMATED in the model needs to be initialised.
inits <- function() list(a=runif(1,35,40), b=runif(1,1,3), c=runif(1,-5,5), 
                         sigma.o=1, sigma.p=1) #small tau makes variance big
inits

function () 
list(a = runif(1, 35, 40), b = runif(1, 1, 3), c = runif(1, -5, 
    5), sigma.o = 1, sigma.p = 1)

# Specify data to be used by your JAGS program (giving dataframe and variables); must be a list. 
data=list(
    n=nrow(hemlock),
    x=as.numeric(hemlock$Light),
    y=as.numeric(hemlock$Observed.growth.rate )
    )
data

$n
[1] 77

$x
 [1] 98.5 98.5 98.5 23.5 11.8 11.8 19.4 19.4 98.7 98.7 98.7 10.9  6.5  6.5 11.1
[16] 50.4 57.5 56.3 53.6 53.5 57.0 67.8 70.9 75.1 22.6 64.3 12.8 29.0 39.6 29.6
[31] 28.7 71.1 98.0 98.0 98.0  5.6  6.9 39.1 28.8 38.2 43.4 38.2 56.4 53.5 29.4
[46] 39.1 18.7 18.0 25.4 57.1 25.9 15.5 14.6 42.9 37.1 66.4 39.3 14.9 14.9 27.4
[61] 19.7 22.9 98.3 98.3 98.3  7.7  7.7  8.8  8.8  8.8 41.3 43.4 40.5 80.5 82.8
[76] 83.5 85.9

$y
 [1] 18.2 26.0 25.7  6.2 10.0 11.2  7.3 15.7 44.2 39.7 30.0  4.8  4.2  5.3  4.7
[16] 28.2 22.7 36.3 37.3 35.3 37.3 33.8 32.0 34.2 12.5 26.2 11.5  7.5 25.0 14.3
[31] 17.2 28.3 12.7 30.3 25.0  2.7  3.5 19.2 30.7 32.2 46.7 24.8 31.7 26.8 29.7
[46] 30.8 28.7 20.2  4.7 25.5  6.8 10.5  4.2 26.3 31.3 13.5 14.3 16.8  8.2 11.7
[61] 23.2 21.2 37.0 34.2 29.2  5.3  6.0  6.5  4.0  6.7 16.8 17.0 23.2 31.7 31.2
[76] 23.5 27.3

# Parameters monitored (i.e., for which estimates are saved)
# params <- c("a","b","c","sigma.o","sigma.p")

# call to JAGS
# Set run conditions: number of iterations for adaptation & runs, number of chains, etc.
n.adapt=500 #plays with algorithm for first 500 samples
n.update = 1000 #burnin
n.iter = 5000 #no. of iterations
jm = jags.model(file = "jags_model.txt",data=data,inits,n.chains=3,n.adapt=n.adapt)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 77
   Unobserved stochastic nodes: 82
   Total graph size: 484

Initializing model

# the jags file name needs to be in the same working directory, or set the path name
# Burnin the chain
update(jm, n.iter=n.update)

The next block of code creates the chains and stores them as an MCMC list. We give the model object name created from when we used the jags.model function, and then variable.names tells JAGS which variables to “monitor” (i.e. the variables for which we want posterior distributions). Then we specify how many iterations we want and to keep - n.thin = 10 means we keep every 10th iteration (which we may choose to do to reduce autocorrelation in our chain).

From summary(objectname) we can get summary statistics for the MCMC chains. The full output of the coda object is a list of matrices, where each matrix is the output of the chains for each parameter (the columns). We can also use the object for plotting (see further below).

# generate coda object for exploring parameters and deviance.
zm <- coda.samples(jm,variable.names=c("a", "b", "c","sigma.o","sigma.p"),n.iter=n.iter, n.thin=10)
summary(zm)

Iterations = 1501:6500
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean     SD Naive SE Time-series SE
a       40.474 5.0478 0.041216        0.24496
b        1.611 0.4587 0.003745        0.02231
c        3.274 2.7177 0.022190        0.12563
sigma.o  5.068 0.5205 0.004250        0.01459
sigma.p  5.397 1.0998 0.008980        0.04462

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
a       32.4815 37.027 39.766 43.228 52.205
b        0.8964  1.288  1.550  1.864  2.675
c       -3.5766  1.861  3.784  5.230  7.069
sigma.o  4.0592  4.711  5.064  5.420  6.113
sigma.p  2.8989  4.773  5.451  6.107  7.413

We also create a jags object, which in brief, are more handy for posterior distributions in matrix form, e.g. posterior distributions of predictions (mu).

# generate JAGS object
zj <- jags.samples(jm,variable.names=c("a", "b","c", "mu2"), n.iter=n.iter, n.thin=10)
zj

$a
mcarray:
[1] 40.79817

Marginalizing over: iteration(5000),chain(3) 

$b
mcarray:
[1] 1.590138

Marginalizing over: iteration(5000),chain(3) 

$c
mcarray:
[1] 3.184362

Marginalizing over: iteration(5000),chain(3) 

$mu2
mcarray:
 [1] 31.286282 31.286282 31.286282 17.304374  9.597032  9.597032 15.119534
 [8] 15.119534 31.301016 31.301016 31.301016  8.759527  3.783739  3.783739
[15]  8.950011 25.531897 26.797382 26.598943 26.130441 26.112473 26.715400
[22] 28.291508 28.679147 29.166648 16.859771 27.822001 10.471764 19.679252
[29] 23.072231 19.907830 19.562955 28.703311 31.249244 31.249244 31.249244
[36]  2.530846  4.310112 22.938046 19.601872 22.691062 24.027002 22.691062
[43] 26.615702 26.112473 19.832224 22.938046 14.701401 14.268111 18.187105
[50] 26.731876 18.407722 12.582421 11.915939 23.907556 22.379268 28.107960
[57] 22.991975 12.141946 12.141946 19.042807 15.294304 17.009975 31.271502
[64] 31.271502 31.271502  5.312057  5.312057  6.590110  6.590110  6.590110
[71] 23.513078 24.027002 23.308515 29.737407 29.963480 30.030422 30.253656

Marginalizing over: iteration(5000),chain(3) 

hat <- summary(zj$mu,quantile,c(0.025,.5,.975))$stat # gives you the median, upper and lower quaniles for mu
hat

          [,1]     [,2]     [,3]     [,4]      [,5]      [,6]     [,7]     [,8]
2.5%  28.16310 28.16310 28.16310 14.81170  6.871919  6.871919 12.63651 12.63651
50%   31.26605 31.26605 31.26605 17.31038  9.610413  9.610413 15.10884 15.10884
97.5% 34.46870 34.46870 34.46870 19.78477 12.328297 12.328297 17.64568 17.64568
          [,9]    [,10]    [,11]     [,12]     [,13]     [,14]     [,15]
2.5%  28.17306 28.17306 28.17306  5.848601 -1.219819 -1.219819  6.091375
50%   31.28137 31.28137 31.28137  8.777812  3.821480  3.821480  8.967320
97.5% 34.49373 34.49373 34.49373 11.656859  8.451086  8.451086 11.807903
         [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]    [,23]
2.5%  23.43272 24.71524 24.52814 24.06027 24.04274 24.63170 26.06238 26.37939
50%   25.55012 26.81406 26.61627 26.14846 26.12884 26.73264 28.30523 28.69390
97.5% 27.54850 28.84404 28.63435 28.15158 28.13371 28.75466 30.51049 30.97761
         [,24]    [,25]    [,26]    [,27]    [,28]    [,29]    [,30]    [,31]
2.5%  26.73317 14.37433 25.66329  7.85258 17.19391 20.79510 17.43397 17.07331
50%   29.17432 16.86426 27.84501 10.48702 19.69509 23.09768 19.92676 19.57782
97.5% 31.57338 19.35230 29.98199 13.09658 22.05515 25.20403 22.26963 21.94344
         [,32]    [,33]    [,34]    [,35]     [,36]      [,37]    [,38]
2.5%  26.39933 28.13927 28.13927 28.13927 -3.363288 -0.3774484 20.64789
50%   28.71887 31.22987 31.22987 31.22987  2.615612  4.3370245 22.96652
97.5% 31.00717 34.41375 34.41375 34.41375  7.840976  8.7492844 25.08081
         [,39]    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
2.5%  17.11310 20.38931 21.83537 20.38931 24.54105 24.04274 17.35365 20.64789
50%   19.61597 22.72152 24.04470 22.72152 26.63207 26.12884 19.84795 22.96652
97.5% 21.98178 24.84657 26.10332 24.84657 28.65251 28.13371 22.20003 25.08081
         [,47]    [,48]    [,49]    [,50]    [,51]    [,52]     [,53]    [,54]
2.5%  12.21241 11.77018 15.69009 24.65083 15.92234 10.10420  9.421289 21.70142
50%   14.69251 14.25556 18.19641 26.74966 18.41652 12.58093 11.917391 23.92429
97.5% 17.21989 16.79088 20.62766 28.77447 20.84061 15.08573 14.433005 25.98919
         [,55]    [,56]    [,57]     [,58]     [,59]    [,60]    [,61]    [,62]
2.5%  20.04902 25.91214 20.70754  9.661187  9.661187 16.54674 12.81301 14.52085
50%   22.40837 28.12446 23.01744 12.144339 12.144339 19.05447 15.28447 17.01677
97.5% 24.55592 30.30616 25.12987 14.644295 14.644295 21.44888 17.81400 19.49381
         [,63]    [,64]    [,65]    [,66]    [,67]     [,68]     [,69]
2.5%  28.15317 28.15317 28.15317 1.205093 1.205093  3.048045  3.048045
50%   31.25051 31.25051 31.25051 5.318754 5.318754  6.594846  6.594846
97.5% 34.44678 34.44678 34.44678 9.320117 9.320117 10.108471 10.108471
          [,70]    [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%   3.048045 21.27688 21.83537 21.05403 27.13823 27.29402 27.34668 27.49788
50%    6.594846 23.53443 24.04470 23.33363 29.73603 29.96686 30.03301 30.25672
97.5% 10.108471 25.61858 26.10332 25.42670 32.29183 32.57241 32.66880 32.98491

# install.packages("lattice")
library(lattice) # this library gives us some additional plotting options
# Plot parameters
# You can set ask to TRUE
plot(zm, ask = FALSE)

plot(zm[, c("a", "b")], ask=FALSE) # subset of parameters

xyplot(zm,ask = FALSE)

densityplot(zm, ask = FALSE)

From these objects we can now plot our predictions along with the data.

# Plot predictions
bu <- summary(zj$mu2, quantile, c(.025,.5,.976))$stat
pl <- cbind(hemlock$Light,t(bu))
dat <- pl[order(pl[,1]),]

plot(hemlock$Light, hemlock$Observed.growth.rate, xlab="Light", ylab="Growth Rate", pch=16, col="blue")
lines(dat[,1], dat[,2],lty="dashed")
lines(dat[,1], dat[,3])
lines(dat[,1], dat[,4],lty="dashed")

We have not covered diagnostics here, but the code below gives some examples of what is possible using the ‘coda’ object you created.

# Convergence diagnostics
rejectionRate(zm) # sampling conjugate

      a       b       c sigma.o sigma.p 
      0       0       0       0       0 

gelman.diag(zm) # var in chains, stable=1

Potential scale reduction factors:

        Point est. Upper C.I.
a             1.01       1.03
b             1.01       1.01
c             1.01       1.01
sigma.o       1.00       1.00
sigma.p       1.00       1.01

Multivariate psrf

1.01

heidel.diag(zm) # requires convergence

[[1]]
                                      
        Stationarity start     p-value
        test         iteration        
a       passed       1         0.209  
b       passed       1         0.511  
c       passed       1         0.533  
sigma.o passed       1         0.121  
sigma.p passed       1         0.150  
                                 
        Halfwidth Mean  Halfwidth
        test                     
a       passed    40.53 0.8203   
b       passed     1.62 0.0848   
c       failed     3.24 0.4004   
sigma.o passed     5.06 0.0449   
sigma.p passed     5.46 0.1265   

[[2]]
                                      
        Stationarity start     p-value
        test         iteration        
a       passed       1         0.261  
b       passed       1         0.126  
c       passed       1         0.105  
sigma.o passed       1         0.640  
sigma.p passed       1         0.519  
                                 
        Halfwidth Mean  Halfwidth
        test                     
a       passed    40.61 0.9358   
b       passed     1.60 0.0729   
c       failed     3.21 0.5140   
sigma.o passed     5.08 0.0439   
sigma.p passed     5.35 0.1403   

[[3]]
                                      
        Stationarity start     p-value
        test         iteration        
a       passed       1         0.399  
b       passed       1         0.438  
c       passed       1         0.554  
sigma.o passed       1         0.172  
sigma.p passed       1         0.285  
                                 
        Halfwidth Mean  Halfwidth
        test                     
a       passed    40.28 0.7252   
b       passed     1.61 0.0687   
c       failed     3.37 0.3481   
sigma.o passed     5.06 0.0584   
sigma.p passed     5.38 0.1821   

raftery.diag(zm) # how many iter you need for convergence

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                               
         Burn-in  Total Lower bound  Dependence
         (M)      (N)   (Nmin)       factor (I)
 a       22       23246 3746          6.21     
 b       48       59565 3746         15.90     
 c       30       30584 3746          8.16     
 sigma.o 7        7968  3746          2.13     
 sigma.p 29       31979 3746          8.54     


[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                                
         Burn-in  Total  Lower bound  Dependence
         (M)      (N)    (Nmin)       factor (I)
 a       13       14172  3746          3.78     
 b       64       66956  3746         17.90     
 c       20       21738  3746          5.80     
 sigma.o 10       10758  3746          2.87     
 sigma.p 116      129088 3746         34.50     


[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                                
         Burn-in  Total  Lower bound  Dependence
         (M)      (N)    (Nmin)       factor (I)
 a       14       15608  3746          4.17     
 b       40       43992  3746         11.70     
 c       30       31743  3746          8.47     
 sigma.o 12       13718  3746          3.66     
 sigma.p 96       101660 3746         27.10     

load.module("dic")

module dic loaded

dic.ex <- dic.samples(jm,n.iter, type="pD")

# autocorrelation with chains
autocorr.plot(zm)

crosscorr(zm)

                  a            b            c      sigma.o      sigma.p
a        1.00000000 -0.794678833 -0.570426234  0.014473938  0.024153685
b       -0.79467883  1.000000000  0.679637251 -0.007489154  0.005789236
c       -0.57042623  0.679637251  1.000000000  0.001083869 -0.065441576
sigma.o  0.01447394 -0.007489154  0.001083869  1.000000000 -0.546326848
sigma.p  0.02415369  0.005789236 -0.065441576 -0.546326848  1.000000000

crosscorr.plot(zm)

Multi-site extension to the model

Now we will explore adding another level to the model - the nature of the relationship between light and plant growth rate can vary according to site, of which we have three in our new dataset. We will model that the maximum growth rate (the asymptote or parameter a) can vary according to site such that this parameter is a random variable with a mean and variance. Lets explore our data briefly again.

jags_model <- "
# JAGS model
########################################################################################
model{
  ### likelihood
  # Data model
    for(s in 1:3){
      for(i in 1:n){
        y[i,s] ~ dnorm(mu[i,s], tau.o) #Note JAGS uses precision, not variance
      }  
    }
  # process model
    for(s in 1:3){
      for(i in 1:n){
        mu[i,s] ~ dnorm(mu2[i,s], tau.p) #Note JAGS uses precision, not variance
    
        mu2[i,s] <- a[s] * (x[i,s]-c) / ((a[s]/b)+(x[i,s]-c))
        }
    }
  
  # priors
  for(s in 1:3){
    a[s] ~ dnorm(mu.a, tau.a)
  }
  mu.a ~ dunif(0,100)
  c ~ dunif(-10,10)
  b ~ dgamma(0.01,0.01)
  
  sigma.o ~ dnorm(5, 1/(0.5*0.5)) ## assume prior knowledge of observation error (5 with SD of 0.5)
  sigma.p ~ dunif(0, 50)
  sigma.a ~ dunif(0, 50)
  
  # derived quantities: convert precisions to standard deviations
  tau.o <- 1/(sigma.o * sigma.o)
  tau.p <- 1/(sigma.p * sigma.p)
  tau.a <- 1/(sigma.a * sigma.a)
  
} # end of model
#########################################################################################
"
writeLines(jags_model, con="jags_model.txt")

# hemlock2 <- read.csv("Hemlock-light-data-hierarchical_2.csv")
# hemlock <- read.csv(url("https://raw.githubusercontent.com/NERC-CEH/beem_data/main/Hemlock-light-data-simple.csv"))
hemlock2 <- read.csv(url("https://raw.githubusercontent.com/NERC-CEH/beem_data/main/Hemlock-light-data-hierarchical_2.csv"))
str(hemlock2) # 231 observations

'data.frame':   231 obs. of  3 variables:
 $ Light               : num  5.6 6.5 6.5 6.9 7.7 7.7 8.8 8.8 8.8 10.9 ...
 $ Observed.growth.rate: num  2.7 4.2 5.3 3.5 5.3 6 6.5 4 6.7 4.8 ...
 $ Site                : int  1 1 1 1 1 1 1 1 1 1 ...

head(hemlock2)

  Light Observed.growth.rate Site
1   5.6                  2.7    1
2   6.5                  4.2    1
3   6.5                  5.3    1
4   6.9                  3.5    1
5   7.7                  5.3    1
6   7.7                  6.0    1

par(mfrow=c(1,1))
colour <- c("red", "blue", "green")
col.list <- rep(0, length(hemlock2$Site))
col.list[hemlock2$Site=="1"] <- 1 
col.list[hemlock2$Site=="2"] <- 2
col.list[hemlock2$Site=="3"] <- 3
plot(hemlock2$Light, hemlock2$Observed.growth.rate, col=c(colour[col.list]))

Again, we need to set up the conditions for the model - how we want to run the MCMC chains, the data etc. Notice how we have to adapt our initial starting values for the new parameteres in the model - open up the new multi-level file.

Note: when you have two product symbols in your conditional distribution with different indices, you have nested for loops.

# Initialize the MCMC chains. This needs to be a list.
inits2 <- function() list(a=rep(30,3), b=runif(1,1,3), c=runif(1,-5,5), 
                         mu.a=30, sigma.o=1, sigma.p=1, sigma.a=1) #small tau makes variance big
inits2

function () 
list(a = rep(30, 3), b = runif(1, 1, 3), c = runif(1, -5, 5), 
    mu.a = 30, sigma.o = 1, sigma.p = 1, sigma.a = 1)

## create dataframe for y data (observed growth rate)

y.dat <- array(NA, c(77,3))
y.dat[1:77,1] <- hemlock2$Observed.growth.rate[which(hemlock2$Site==1)]
y.dat[1:77,2] <- hemlock2$Observed.growth.rate[which(hemlock2$Site==2)]
y.dat[1:77,3] <- hemlock2$Observed.growth.rate[which(hemlock2$Site==3)]
head(y.dat)

     [,1] [,2] [,3]
[1,]  2.7  4.7  2.6
[2,]  4.2  2.7  5.5
[3,]  5.3  5.9  5.2
[4,]  3.5  2.3  4.3
[5,]  5.3  4.7  4.2
[6,]  6.0  7.0  6.6

## create dataframe for x data (Light)
x.dat <- array(NA, c(77,3))
x.dat[1:77,1] <- hemlock2$Light[which(hemlock2$Site==1)]
x.dat[1:77,2] <- hemlock2$Light[which(hemlock2$Site==2)]
x.dat[1:77,3] <- hemlock2$Light[which(hemlock2$Site==3)]
head(x.dat)

     [,1] [,2] [,3]
[1,]  5.6  3.9  5.5
[2,]  6.5  4.8  5.7
[3,]  6.5  6.7  5.9
[4,]  6.9  7.1  6.5
[5,]  7.7  7.6  7.8
[6,]  7.7  8.3  7.8

# Specify data; must be a list. 
data=list(
  n=77,
  x=x.dat,
  y=y.dat
)
data

$n
[1] 77

$x
      [,1] [,2]  [,3]
 [1,]  5.6  3.9   5.5
 [2,]  6.5  4.8   5.7
 [3,]  6.5  6.7   5.9
 [4,]  6.9  7.1   6.5
 [5,]  7.7  7.6   7.8
 [6,]  7.7  8.3   7.8
 [7,]  8.8  8.7   8.0
 [8,]  8.8  9.2   8.3
 [9,]  8.8  9.3   9.2
[10,] 10.9  9.3  10.7
[11,] 11.1  9.4  11.1
[12,] 11.8 10.4  11.6
[13,] 11.8 11.5  12.2
[14,] 12.8 11.5  13.5
[15,] 14.6 13.1  13.8
[16,] 14.9 14.0  15.1
[17,] 14.9 14.9  15.6
[18,] 15.5 16.0  17.4
[19,] 18.0 17.5  17.5
[20,] 18.7 18.1  17.9
[21,] 19.4 18.4  18.6
[22,] 19.4 20.0  18.9
[23,] 19.7 21.0  20.6
[24,] 22.6 21.1  21.0
[25,] 22.9 21.5  22.4
[26,] 23.5 22.9  23.9
[27,] 25.4 24.0  24.3
[28,] 25.9 24.3  25.1
[29,] 27.4 27.1  28.3
[30,] 28.7 27.3  28.6
[31,] 28.8 28.9  28.9
[32,] 29.0 29.0  29.1
[33,] 29.4 30.2  30.5
[34,] 29.6 31.2  30.7
[35,] 37.1 36.2  36.4
[36,] 38.2 38.2  37.4
[37,] 38.2 38.4  37.8
[38,] 39.1 38.8  38.2
[39,] 39.1 39.7  39.8
[40,] 39.3 40.1  39.8
[41,] 39.6 40.5  40.1
[42,] 40.5 40.7  40.8
[43,] 41.3 41.0  41.4
[44,] 42.9 41.8  41.8
[45,] 43.4 42.8  42.7
[46,] 43.4 44.2  43.7
[47,] 50.4 49.4  51.5
[48,] 53.5 52.0  51.8
[49,] 53.5 53.6  52.9
[50,] 53.6 54.6  55.1
[51,] 56.3 55.9  55.5
[52,] 56.4 57.3  56.8
[53,] 57.0 57.9  57.2
[54,] 57.1 58.0  57.5
[55,] 57.5 58.7  58.3
[56,] 64.3 64.9  62.5
[57,] 66.4 66.2  65.1
[58,] 67.8 66.6  68.5
[59,] 70.9 69.1  69.7
[60,] 71.1 72.2  70.2
[61,] 75.1 76.1  74.3
[62,] 80.5 80.8  81.2
[63,] 82.8 82.0  82.2
[64,] 83.5 84.1  84.7
[65,] 85.9 85.0  86.6
[66,] 98.0 96.8  96.2
[67,] 98.0 97.0  96.8
[68,] 98.0 97.1  97.2
[69,] 98.3 97.5  97.3
[70,] 98.3 98.2  97.6
[71,] 98.3 98.5  98.2
[72,] 98.5 99.4  98.5
[73,] 98.5 99.4  98.8
[74,] 98.5 99.5  98.9
[75,] 98.7 99.6  99.7
[76,] 98.7 99.7 100.0
[77,] 98.7 99.9 100.3

$y
      [,1]      [,2]     [,3]
 [1,]  2.7  4.700000  2.60000
 [2,]  4.2  2.700000  5.50000
 [3,]  5.3  5.900000  5.20000
 [4,]  3.5  2.300000  4.30000
 [5,]  5.3  4.700000  4.20000
 [6,]  6.0  7.000000  6.60000
 [7,]  6.5  7.000000  4.30000
 [8,]  4.0  2.700000  3.00000
 [9,]  6.7  5.000000  6.50000
[10,]  4.8  5.800000  5.30000
[11,]  4.7  3.200000  3.00000
[12,] 10.0  8.200000 10.60000
[13,] 11.2  9.300000 12.30000
[14,] 11.5  9.700000  9.90000
[15,]  4.2  3.200000  1.40000
[16,] 16.8 15.200000 17.60000
[17,]  8.2  6.700000  7.70000
[18,] 10.5 11.500000  8.60000
[19,] 20.2 19.900000 20.20000
[20,] 28.7 27.700000 28.50000
[21,]  7.3  7.800000  6.50000
[22,] 15.7 15.300000 16.90000
[23,] 23.2 22.800000 22.80000
[24,] 12.5 13.000000 14.80000
[25,] 21.2 22.600000 22.20000
[26,]  6.2  7.700000  9.60000
[27,]  4.7  3.900000  5.50000
[28,]  6.8  5.400000 10.60000
[29,] 11.7 12.300000  8.70000
[30,] 17.2 15.500000 18.10000
[31,] 30.7 32.500000 33.70000
[32,]  7.5  6.700000  5.00000
[33,] 29.7 31.300000 29.70000
[34,] 14.3 15.700000 12.50000
[35,] 31.3 47.136898 63.89346
[36,] 32.2 28.957552 53.82420
[37,] 24.8 33.877043 85.72795
[38,] 19.2 32.241437 69.38528
[39,] 30.8 23.643437 35.94980
[40,] 14.3 39.105755 40.44215
[41,] 25.0 43.813701 37.19800
[42,] 23.2 32.523262 46.49454
[43,] 16.8 21.814450 49.02929
[44,] 26.3 30.259437 50.90613
[45,] 46.7 18.943678 66.51076
[46,] 17.0 45.630525 52.30740
[47,] 28.2 20.099826 90.76589
[48,] 35.3 18.316243 87.56340
[49,] 26.8 15.446236 39.93826
[50,] 37.3  8.010457 50.31908
[51,] 36.3 31.162236 67.13457
[52,] 31.7 31.934466 59.57232
[53,] 37.3 16.407239 50.69410
[54,] 25.5 19.764969 55.26655
[55,] 22.7 41.048627 56.47245
[56,] 26.2 34.287661 67.75295
[57,] 13.5 31.308202 76.56862
[58,] 33.8 16.114231 40.12719
[59,] 32.0 32.443813 48.78845
[60,] 28.3 32.601453 41.05356
[61,] 34.2 37.722595 50.32160
[62,] 31.7 28.820051 57.95963
[63,] 31.2 17.373145 56.57387
[64,] 23.5 31.298986 54.23894
[65,] 27.3 28.693491 67.87955
[66,] 12.7 39.714315 26.83006
[67,] 30.3 39.255526 39.53883
[68,] 25.0 36.226181 52.91146
[69,] 37.0 15.376946 60.52105
[70,] 34.2 44.137594 56.51756
[71,] 29.2 29.533540 71.09812
[72,] 18.2 37.031727 71.52053
[73,] 26.0 21.899740 63.90977
[74,] 25.7 21.608192 69.75787
[75,] 44.2 24.773272 43.06696
[76,] 39.7 45.428707 31.40532
[77,] 30.0 21.766299 57.56964

# Parameters monitored (i.e., for which estimates are saved)
# params <- c("a","b","c","sigma.o","sigma.p", "sigma.a", "mu.a")

Then we call to JAGS using our model file and the rjags function jags.model.

### fit the model in JAGS
n.adapt = 500
n.update = 1000
n.iter = 5000
jm2 = jags.model(file = "jags_model.txt", data = data, inits2, n.chains = 3, n.adapt=n.adapt)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 231
   Unobserved stochastic nodes: 240
   Total graph size: 1528

Initializing model

# Burn in the chain
update(jm2, n.iter=n.update)

# Generate coda object
zm2 <- coda.samples(jm2,variable.names=c("a","b","c","sigma.o","sigma.p","sigma.a", "mu.a"), n.iter=n.iter,thin=1)
summary(zm)

Iterations = 1501:6500
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean     SD Naive SE Time-series SE
a       40.474 5.0478 0.041216        0.24496
b        1.611 0.4587 0.003745        0.02231
c        3.274 2.7177 0.022190        0.12563
sigma.o  5.068 0.5205 0.004250        0.01459
sigma.p  5.397 1.0998 0.008980        0.04462

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
a       32.4815 37.027 39.766 43.228 52.205
b        0.8964  1.288  1.550  1.864  2.675
c       -3.5766  1.861  3.784  5.230  7.069
sigma.o  4.0592  4.711  5.064  5.420  6.113
sigma.p  2.8989  4.773  5.451  6.107  7.413

# Generate jags object
zj2 <- jags.samples(jm2,variable.names=c("a", "b","c", "mu2"), n.iter=n.iter, n.thin=10)
zj2

$a
mcarray:
[1] 36.09709 37.49914 88.49663

Marginalizing over: iteration(5000),chain(3) 

$b
mcarray:
[1] 2.225225

Marginalizing over: iteration(5000),chain(3) 

$c
mcarray:
[1] 5.804087

Marginalizing over: iteration(5000),chain(3) 

$mu2
mcarray:
            [,1]      [,2]        [,3]
 [1,] -0.7301526 -5.190311 -0.86930587
 [2,]  1.2627876 -2.697695 -0.41548493
 [3,]  1.2627876  1.682110  0.03355998
 [4,]  2.0784769  2.481124  1.35280992
 [5,]  3.5967801  3.429512  4.07551373
 [6,]  3.5967801  4.671572  4.07551373
 [7,]  5.4698959  5.340326  4.47880431
 [8,]  5.4698959  6.137868  5.07630447
 [9,]  5.4698959  6.292523  6.81713195
[10,]  8.4930770  6.292523  9.55739242
[11,]  8.7494533  6.445619 10.25636585
[12,]  9.6096660  7.895886 11.11232766
[13,]  9.6096660  9.339056 12.11433705
[14,] 10.7469348  9.339056 14.19627558
[15,] 12.5616974 11.199114 14.66025112
[16,] 12.8387804 12.139404 16.60386126
[17,] 12.8387804 13.013686 17.32371957
[18,] 13.3735212 14.002482 19.79717310
[19,] 15.3568881 15.227153 19.92943307
[20,] 15.8509044 15.681537 20.45329761
[21,] 16.3216571 15.901756 21.35060751
[22,] 16.3216571 17.004165 21.72779782
[23,] 16.5166945 17.637170 23.78594532
[24,] 18.2191919 17.698303 24.25151528
[25,] 18.3784027 17.939076 25.82849888
[26,] 18.6883875 18.737109 27.43279207
[27,] 19.6016344 19.319457 27.84656148
[28,] 19.8261249 19.472000 28.65719358
[29,] 20.4641046 20.780341 31.69057073
[30,] 20.9776023 20.866507 31.95899575
[31,] 21.0156985 21.524904 32.22485419
[32,] 21.0913137 21.564315 32.40068569
[33,] 21.2402765 22.022411 33.60092366
[34,] 21.3136468 22.384379 33.76814087
[35,] 23.6214891 23.965635 38.13457374
[36,] 23.9004602 24.508785 38.82885101
[37,] 23.9004602 24.560700 39.10116318
[38,] 24.1195545 24.663289 39.37046797
[39,] 24.1195545 24.888243 40.41859952
[40,] 24.1671726 24.985700 40.41859952
[41,] 24.2378909 25.081656 40.61010506
[42,] 24.4450819 25.129081 41.05100330
[43,] 24.6232419 25.199541 41.42242990
[44,] 24.9636395 25.383552 41.66679648
[45,] 25.0658939 25.605933 42.20734864
[46,] 25.0658939 25.903843 42.79335225
[47,] 26.3221689 26.890426 46.89649443
[48,] 26.7912477 27.323144 47.03951153
[49,] 26.7912477 27.572266 47.55552505
[50,] 26.8056287 27.721849 48.54963195
[51,] 27.1777480 27.909679 48.72517117
[52,] 27.1909572 28.104030 49.28508286
[53,] 27.2693957 28.184926 49.45417967
[54,] 27.2823343 28.198273 49.58004055
[55,] 27.3337108 28.290639 49.91170485
[56,] 28.1235686 29.034974 51.56351897
[57,] 28.3394445 29.176006 52.51654675
[58,] 28.4769707 29.218458 53.69055581
[59,] 28.7647666 29.474257 54.08671277
[60,] 28.7825881 29.770271 54.24911240
[61,] 29.1216409 30.112983 55.52479535
[62,] 29.5323721 30.487602 57.46872857
[63,] 29.6929647 30.577229 57.73150484
[64,] 29.7402906 30.728675 58.36928214
[65,] 29.8973653 30.791568 58.83649859
[66,] 30.5840502 31.519705 60.99083801
[67,] 30.5840502 31.530685 61.11509746
[68,] 30.5840502 31.536160 61.19731270
[69,] 30.5991463 31.557959 61.21778910
[70,] 30.5991463 31.595726 61.27903369
[71,] 30.5991463 31.611765 61.40069820
[72,] 30.6091645 31.659363 61.46112185
[73,] 30.6091645 31.659363 61.52127553
[74,] 30.6091645 31.664605 61.54126707
[75,] 30.6191462 31.669836 61.70013480
[76,] 30.6191462 31.675059 61.75922658
[77,] 30.6191462 31.685475 61.81805728

Marginalizing over: iteration(5000),chain(3) 

# summary(zj2$mu.a, quantile, c(0.025, 0.5, 0.975))$stat

# Plot parameters
# library(lattice) required for these plots
xyplot(zm2,ask = TRUE)

densityplot(zm2,ask = TRUE)

summary(zm2)

Iterations = 1501:6500
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean      SD Naive SE Time-series SE
a[1]    36.094  3.0266 0.024712       0.054364
a[2]    37.443  3.0746 0.025104       0.054247
a[3]    88.678  7.6146 0.062173       0.216031
b        2.220  0.2683 0.002191       0.007823
c        5.773  1.0186 0.008317       0.021164
mu.a    53.800 18.0960 0.147753       0.207976
sigma.a 32.515  9.3075 0.075996       0.113203
sigma.o  5.026  0.5095 0.004160       0.016010
sigma.p  9.954  0.6475 0.005287       0.014207

2. Quantiles for each variable:

          2.5%    25%    50%    75%   97.5%
a[1]    30.454 34.008 35.995 38.048  42.362
a[2]    31.760 35.322 37.342 39.463  43.622
a[3]    75.685 83.318 87.978 93.293 105.241
b        1.731  2.033  2.208  2.394   2.776
c        3.486  5.166  5.877  6.491   7.478
mu.a    16.398 42.120 53.965 65.601  89.785
sigma.a 15.969 25.191 32.287 40.004  48.825
sigma.o  4.009  4.688  5.028  5.367   6.013
sigma.p  8.703  9.521  9.949 10.382  11.244

summary(zm2)$stat

             Mean         SD    Naive SE Time-series SE
a[1]    36.094073  3.0266153 0.024712211    0.054363587
a[2]    37.443439  3.0746134 0.025104113    0.054246817
a[3]    88.678431  7.6146461 0.062173325    0.216030527
b        2.219596  0.2683226 0.002190845    0.007823399
c        5.773428  1.0185932 0.008316778    0.021164220
mu.a    53.800137 18.0959758 0.147753024    0.207976284
sigma.a 32.514631  9.3075300 0.075995664    0.113202531
sigma.o  5.026196  0.5094796 0.004159883    0.016010458
sigma.p  9.954361  0.6474827 0.005286674    0.014206717

# Plot predictions
bu2 <- summary(zj2$mu2,quantile,c(.025,.5,.976))$stat
bu2

, , 1

            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]     [,7]
2.5%  -5.3259874 -2.584404 -2.584404 -1.490844 0.5009855 0.5009855 2.844969
50%   -0.6553475  1.301018  1.301018  2.109491 3.6185951 3.6185951 5.481122
97.6%  3.5486892  4.912091  4.912091  5.498218 6.6096403 6.6096403 8.051871
          [,8]     [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
2.5%  2.844969 2.844969  6.393317  6.688836  7.642421  7.642421  8.867437
50%   5.481122 5.481122  8.496797  8.752090  9.614769  9.614769 10.752528
97.6% 8.051871 8.051871 10.584399 10.806002 11.569362 11.569362 12.608273
         [,15]    [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]
2.5%  10.79297 11.09478 11.09478 11.64445 13.65143 14.14404 14.60443 14.60443
50%   12.55550 12.83090 12.83090 13.36411 15.33984 15.83681 16.30601 16.30601
97.6% 14.34466 14.60976 14.60976 15.13840 17.13370 17.63696 18.11531 18.11531
         [,23]    [,24]    [,25]    [,26]    [,27]    [,28]    [,29]    [,30]
2.5%  14.79780 16.42174 16.56734 16.86558 17.71125 17.91091 18.48999 18.94909
50%   16.49876 18.20618 18.36503 18.67750 19.59097 19.81648 20.45432 20.97244
97.6% 18.32223 20.09407 20.26090 20.58670 21.56345 21.80576 22.48890 23.05446
         [,31]    [,32]    [,33]    [,34]    [,35]    [,36]    [,37]    [,38]
2.5%  18.98280 19.05418 19.18647 19.25047 21.25618 21.49461 21.49461 21.67518
50%   21.01093 21.08540 21.23443 21.30828 23.61582 23.89504 23.89504 24.11554
97.6% 23.09773 23.17902 23.34605 23.42959 26.02701 26.34863 26.34863 26.59911
         [,39]    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
2.5%  21.67518 21.71408 21.77130 21.94462 22.10240 22.38960 22.47992 22.47992
50%   24.11554 24.16427 24.23381 24.44025 24.61682 24.95275 25.05507 25.05507
97.6% 26.59911 26.65587 26.73509 26.98060 27.17323 27.57573 27.69338 27.69338
         [,47]    [,48]    [,49]    [,50]    [,51]    [,52]    [,53]    [,54]
2.5%  23.50542 23.87766 23.87766 23.88981 24.17324 24.18263 24.24459 24.25484
50%   26.31082 26.78254 26.78254 26.79661 27.16825 27.18121 27.25959 27.27260
97.6% 29.18789 29.76047 29.76047 29.78031 30.22304 30.23874 30.33019 30.34665
         [,55]    [,56]    [,57]    [,58]    [,59]    [,60]    [,61]    [,62]
2.5%  24.29552 24.91569 25.08421 25.19134 25.40442 25.41637 25.68281 25.98345
50%   27.32414 28.11247 28.32669 28.46444 28.75459 28.77302 29.11446 29.51627
97.6% 30.40926 31.37968 31.65764 31.82957 32.18924 32.20897 32.62748 33.15175
         [,63]    [,64]    [,65]    [,66]    [,67]    [,68]    [,69]    [,70]
2.5%  26.10731 26.14337 26.25706 26.77420 26.77420 26.77420 26.78720 26.78720
50%   29.67659 29.72435 29.88075 30.55944 30.55944 30.55944 30.57445 30.57445
97.6% 33.35622 33.42143 33.62890 34.53625 34.53625 34.53625 34.55790 34.55790
         [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%  26.78720 26.79582 26.79582 26.79582 26.80403 26.80403 26.80403
50%   30.57445 30.58414 30.58414 30.58414 30.59444 30.59444 30.59444
97.6% 34.55790 34.57070 34.57070 34.57070 34.58340 34.58340 34.58340

, , 2

             [,1]      [,2]      [,3]       [,4]      [,5]     [,6]     [,7]
2.5%  -11.9172354 -8.089491 -2.019042 -0.9780419 0.2567641 1.838114 2.657276
50%    -5.0256550 -2.609294  1.716920  2.5111202 3.4521011 4.686058 5.353604
97.6%   0.6159773  2.246469  5.232968  5.8227982 6.5203141 7.465516 7.985196
          [,8]     [,9]    [,10]    [,11]     [,12]     [,13]     [,14]
2.5%  3.616148 3.795222 3.795222 3.967204  5.653886  7.300380  7.300380
50%   6.149762 6.302098 6.302098 6.452257  7.901697  9.340906  9.340906
97.6% 8.620211 8.746352 8.746352 8.878560 10.105885 11.376996 11.376996
          [,15]    [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]
2.5%   9.303083 10.29873 11.21248 12.24197 13.47685 13.93239 14.15475 15.23714
50%   11.195787 12.13418 13.00737 13.99615 15.22217 15.67212 15.89121 16.99315
97.6% 13.078717 13.96882 14.82657 15.80862 17.02957 17.50422 17.72419 18.84699
         [,23]    [,24]    [,25]    [,26]    [,27]    [,28]    [,29]    [,30]
2.5%  15.85240 15.91397 16.14189 16.89606 17.42753 17.57210 18.76353 18.84270
50%   17.62708 17.68976 17.93096 18.73798 19.32050 19.47168 20.78350 20.86989
97.6% 19.50161 19.56952 19.82219 20.65687 21.26510 21.41963 22.83233 22.92759
         [,31]    [,32]    [,33]    [,34]    [,35]    [,36]    [,37]    [,38]
2.5%  19.44027 19.47613 19.87818 20.19951 21.56472 22.03555 22.07760 22.16272
50%   21.52825 21.56711 22.02692 22.39112 23.97812 24.52519 24.57684 24.68065
97.6% 23.61778 23.66258 24.15327 24.55160 26.30881 26.93167 26.99334 27.10586
         [,39]    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
2.5%  22.35232 22.43442 22.51517 22.55507 22.61504 22.76935 22.95545 23.20528
50%   24.90574 25.00247 25.09758 25.14578 25.21562 25.40165 25.62750 25.92460
97.6% 27.36562 27.48716 27.59627 27.64534 27.72865 27.93204 28.18383 28.53704
         [,47]    [,48]    [,49]    [,50]    [,51]    [,52]    [,53]    [,54]
2.5%  24.01885 24.36535 24.56087 24.67975 24.82987 24.98499 25.05212 25.06408
50%   26.91438 27.34396 27.59291 27.74281 27.93014 28.12277 28.20427 28.21685
97.6% 29.70536 30.22576 30.52772 30.70918 30.93870 31.17901 31.27124 31.28658
         [,55]    [,56]    [,57]    [,58]    [,59]    [,60]    [,61]    [,62]
2.5%  25.14403 25.72893 25.84279 25.87679 26.07807 26.29511 26.55291 26.82743
50%   28.31071 29.04823 29.18512 29.22573 29.47990 29.77610 30.11401 30.48208
97.6% 31.40831 32.33882 32.51060 32.56455 32.89309 33.27352 33.70835 34.19590
         [,63]    [,64]    [,65]    [,66]    [,67]    [,68]    [,69]    [,70]
2.5%  26.89669 27.01469 27.05834 27.59456 27.60166 27.60520 27.62049 27.64857
50%   30.56823 30.71805 30.77885 31.50321 31.51396 31.52002 31.54231 31.57893
97.6% 34.30896 34.50174 34.58751 35.55415 35.56769 35.57444 35.60296 35.65391
         [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%  27.65973 27.69292 27.69292 27.69657 27.70022 27.70417 27.71327
50%   31.59432 31.64222 31.64222 31.64715 31.65219 31.65741 31.66818
97.6% 35.67734 35.74076 35.74076 35.74850 35.75622 35.76393 35.77931

, , 3

            [,1]       [,2]        [,3]      [,4]      [,5]      [,6]     [,7]
2.5%  -5.1745207 -4.6096910 -4.05937872 -2.480284 0.7419107 0.7419107 1.212171
50%   -0.8746799 -0.4255879  0.01413717  1.328667 4.0526772 4.0526772 4.453684
97.6%  3.5948489  3.9424410  4.30798992  5.357614 7.6088818 7.6088818 7.947429
          [,8]     [,9]     [,10]     [,11]     [,12]     [,13]    [,14]
2.5%  1.892229 3.823604  6.759219  7.492383  8.368824  9.400825 11.50596
50%   5.054543 6.796829  9.527199 10.233062 11.084236 12.091058 14.17922
97.6% 8.446176 9.984504 12.472718 13.138604 13.948718 14.914133 16.94440
         [,15]    [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]
2.5%  11.97830 13.90068 14.62514 17.09631 17.22793 17.74835 18.62898 18.99000
50%   14.64341 16.58394 17.30692 19.77301 19.90383 20.42775 21.32799 21.70372
97.6% 17.40246 19.39181 20.12218 22.66091 22.79318 23.32921 24.25610 24.62785
         [,23]    [,24]    [,25]    [,26]    [,27]    [,28]    [,29]    [,30]
2.5%  20.99612 21.45381 22.99024 24.55093 24.96322 25.75628 28.73400 28.99779
50%   23.75884 24.22416 25.80424 27.41142 27.82477 28.63374 31.67167 31.94246
97.6% 26.72152 27.20348 28.81338 30.44319 30.86246 31.68273 34.72542 34.99927
         [,31]    [,32]    [,33]    [,34]    [,35]    [,36]    [,37]    [,38]
2.5%  29.26625 29.43701 30.63315 30.80244 35.16801 35.86520 36.13563 36.40657
50%   32.20887 32.38379 33.58443 33.75211 38.11809 38.81671 39.09071 39.36305
97.6% 35.26258 35.43716 36.62541 36.78871 41.13245 41.80767 42.07617 42.33152
         [,39]    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
2.5%  37.46580 37.46580 37.65564 38.09715 38.47560 38.71459 39.25537 39.85168
50%   40.41241 40.41241 40.60250 41.04112 41.41217 41.65512 42.19319 42.78169
97.6% 43.36626 43.36626 43.55804 44.00178 44.38418 44.62868 45.16214 45.74758
         [,47]    [,48]    [,49]    [,50]    [,51]    [,52]    [,53]    [,54]
2.5%  43.94331 44.08077 44.59166 45.54908 45.71125 46.26087 46.41595 46.53479
50%   46.89511 47.03709 47.55416 48.54163 48.71740 49.27346 49.44069 49.56686
97.6% 49.88328 50.02880 50.56868 51.59934 51.78496 52.36717 52.54706 52.68665
         [,55]    [,56]    [,57]    [,58]    [,59]    [,60]    [,61]    [,62]
2.5%  46.85180 48.39456 49.25593 50.32252 50.68667 50.82786 51.98064 53.69766
50%   49.90219 51.54772 52.50302 53.67498 54.07506 54.23462 55.51533 57.44978
97.6% 53.02850 54.77216 55.81053 57.10575 57.54758 57.72668 59.16562 61.35716
         [,63]    [,64]    [,65]    [,66]    [,67]    [,68]    [,69]    [,70]
2.5%  53.92685 54.46134 54.84235 56.60461 56.70899 56.77248 56.78679 56.83617
50%   57.71386 58.35408 58.81148 60.95984 61.08348 61.16384 61.18425 61.24657
97.6% 61.66153 62.41339 62.96435 65.54653 65.69386 65.78842 65.81194 65.88393
         [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%  56.93657 56.98048 57.02835 57.04424 57.17067 57.21786 57.26826
50%   61.36786 61.42938 61.48767 61.50794 61.66535 61.72474 61.78242
97.6% 66.03546 66.10652 66.18340 66.20754 66.39869 66.47546 66.55084

Additional resources

We have not used this package here, but do take a look at the functionality offered by the package ‘jagshelper’ - https://cran.r-project.org/web/packages/jagshelper/vignettes/jagshelper-vignette.html

install.packages("jagshelper")
library(jagshelper)

Plummer (2017) JAGS user manual - https://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf